function [thr, fa, fr] = slgetroc(thrs, fars, frrs, item, itempara)
%SLGETROC Computes some point from ROC Curve
%
% [ Syntax ]
%   - [thr, fa, fr] = slgetroc(thrs, fars, frrs, item, itempara)
%
% [ Arguments ]
%   - thrs:         the sampled threshold values
%   - fars:         the false accept rates at the sampled thresholds
%   - frrs:         the false reject rates at the sampled thresholds
%   - item:         the items to be evaluated
%   - itempara:     the extra parameter for the item
%
%   - thr:          the threshold at which the specified point is found
%   - fa:           the false accept rate at the selected point
%   - fr:           the false reject rate at the selected point
%   
% [ Description ]
%   [thr, fa, fr] = slgetroc(thrs, fars, frrs, item, itempara) finds 
%   a point that satisfies a specified constraint from the ROC curves 
%   given by thrs, fars, and frrs. 
%
%   The requirement on the point is specified by item and itempara.
%   \{:
%      - ratio:     solves the point where fr / fa == itempara
%      - fixfa:     solves the point where fa == itempara
%      - fixfr:     solves the point where fr == itempara
%      - fixth:     solves the point where threshold == itempara 
%   \:}
%   The item can also be a function handle with the syntax
%       v = f(thr, fa, fr)
%   which produces a monotonical vector of values, from which the function
%   finds the point v == itempara.
%
%   If itempara is a scalar, then all outputs are scalars. If itempara is
%   an array, then all outputs are arrays of the same size as itempara.
%   
% [ Remarks ]
%   - To increase accuracy, interpolation technique is used.
%
% [ History ]
%   - Created by Dahua Lin on Jun 10th, 2005
%   - Modified by Dahua Lin on May 1st, 2006
%   - Modified by Dahua Lin on Jul 19, 2007
%       - reimplement the interpolation
%       - refactor the computation
%       - support user-supplied item
%       - support multiple item parameter
% 

%% Parse and verify input arguments

error(nargchk(5, 5, nargin));

assert(isnumeric(thrs) && isvector(thrs) && length(thrs) >= 2, ...
    'sltoolbox:slgetroc:invalidarg', 'thrs should be a numeric vector with length not less than 2.');
assert(isnumeric(fars) && isvector(fars) && length(fars) >= 2, ...
    'sltoolbox:slgetroc:invalidarg', 'fars should be a numeric vector with length not less than 2.');
assert(isnumeric(frrs) && isvector(frrs) && length(frrs) >= 2, ...
    'sltoolbox:slgetroc:invalidarg', 'frrs should be a numeric vector with length not less than 2.');

assert(isequal(size(thrs), size(fars), size(frrs)), 'sltoolbox:slgetroc:invalidarg', ...
    'The sizes of thrs, fars, and frrs should be the same.');

assert(ischar(item) || isa(item, 'function_handle'), 'sltoolbox:slgetroc:invalidarg', ...
    'item should be a string or a function handle.');

assert(isnumeric(itempara), 'sltoolbox:slgetroc:invalidarg', ...
    'itempara should be a numeric array.');


%% main

Y0 = [thrs(:) fars(:) frrs(:)];

switch item
    case 'ratio'
        x0 = Y0(:,3) ./ Y0(:,2);                
    case 'fixfa'
        x0 = Y0(:,2);        
    case 'fixfr'
        x0 = Y0(:,3);        
    case 'fixth'
        x0 = Y0(:,1);        
    otherwise
        error('sltoolbox:invalidarg', ...
            'Invalid item %s for slgetroc', item);
end

Yq = my_interp(x0, Y0, itempara(:));

thr = Yq(:, 1);
fa = Yq(:, 2);
fr = Yq(:, 3);

out_size = size(itempara);
if ~isequal(size(thr), out_size)
    thr = reshape(thr, out_size);
    fa = reshape(fa, out_size);
    fr = reshape(fr, out_size);
end


%% Specific interpolation implementation

function Yq = my_interp(x0, Y0, xq)

% reduce the references

[sps, eps, xr] = sllabelrgns(x0, 'start','end', 'val');

if length(xr) == length(x0)
    Yr = Y0;
else
    Yr = Y0(sps, :);
    coll = (eps > sps);
    cY = arrayfun(@(sp, ep) sum(Y0(sp:ep, :), 1) / (ep-sp+1), ...
        sps(coll), eps(coll), 'UniformOutput', false);
    Yr(coll, :) = vertcat(cY{:});
end

% confine the query

min_x = min(x0);
max_x = max(x0);

xq(xq < min_x) = min_x;
xq(xq > max_x) = max_x;

% do interpolation

Yq = interp1(xr, Yr, xq);

